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Abstract 

Q We develop a numerical scheme for solving time-domain Maxwell's equation. The method 

is motivated by CIP method which uses function values and its derivatives as unknown 
variables. The proposed scheme is developed by using the Poisson formula for the wave 
equation. It is fully explicit space and time integration method with higher order accuracy 
and CFL number being one. The bi-cubic interpolation is used for the solution profile to 
I— I attain the resolution. It preserves sharp profiles very accurately without any smearing and 

distortion due to the exact time integration and high resolution approximation. The stability 
and numerical accuracy are investigated. 

c3 1 Introduction 

s 

'— ' Multi-moment methods for time dependent differential equations aim to increase the accuracy 

of the numerical solution, and to lower the dispersive and dissipative errors in the numerical 
J> solution. The most distinguishing characteristic of the method is that more than one moment 

per grid point or per cell, such as function value, its derivatives and the integral over the cell 
etc, are considered as unknown variables, and they are simultaneously updated by coupling the 

•/^ differential equation and derived differential equations for the derivatives. A Hermite polynomial 

is usually defined on each cell using such quantities to interpolate the numerical solution, and 

^ is used to update the solutions. Such polynomial interpolation, being defined on each cell 

individually, reduces the numerical stencil width. The compactness of the stencil makes it 
feasible to handle the boundary condition and the interface condition numerically. 

• Such schemes were first proposed for one dimensional hyperbolic equations by Van-Leer [9] in 

rS a framework of a finite difference method, and by Takewaki, Nishiguchi and Yabe [15] on the basis 

of the characteristic equation for the first order PDE. The latter method is referred to as CIP 
method. Their works triggered off development of multi moment methods, and it is increasingly 
becoming an active area of research and have even been applied to various equations. There is 
a vast literature developing CIP methods and the following is a partial list of papers: nonlinear 
hyperbolic equations [17J, multi dimensional hyperbolic equations fl6l I18j. a multi-dimensional 
the Maxwell's equations [13j, a numerical simulation for solid, liquid and gas [12l [19l [20] , a new 
mesh system applicable to non-orthogonal coordinate system [21], a variant of CIP method [2]. 
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CIP method uses the cubic interpolation constructed via solution values and its derivative 
at two end points of a cell to approximate the solution in the cell. Explicit time integration 
formulas for both the exact solution and its derivative of a one dimensional transport equation 
with the constant velocity field are coupled with the cubic interpolation to obtain an explicit 
time integration numerical scheme (CIP scheme). In [8] we have developed and analyzed a CIP 
scheme for one dimensional hyperbolic equations with variable and discontinuous coefficients. 

Most of the proposed CIP methods so far, when it is applied to higher dimensional equations, 
relies on the Strang splitting technique: It reduces the equation of higher dimension to a sequence 
of simpler one dimensional equations, repeatedly applies a CIP method for the reduced one 
dimensional equations. 

In this paper we develop a multi-moment method for two dimensional Maxwell's equations, 
which does not employ the directional splitting technique at all. We apply the exact integration 
method in time by the Poisson formula and use the bi-cubic interpolation. We refer [3j [lOj 
[Z!)[I1]; m for numerical integration methods for the Maxwell and acoustic wave equations 
in general. The original second order Yee scheme using the time and space staggered grid is 
developed in [22] and a fourth order time and space variant of Yee scheme is developed in [31 E] • 

Our contributions are as follows. The bi-cubic interpolation is combined with the Poisson 
formula to develop a fully explicit time and space numerical method for Maxwell, acoustic wave 
equations as well as the second order wave equation. We analyze the von-Neumann stability 
of the proposed method and we establish CFL number is one for the proposed method. The 
one step of the method involves updating four moments at each grid point and the symbolic 
formulation and the symmetry of the update is used to reduce operation counts. The method 
offers highly efficient to attain the desired accuracy due to the relaxed CFL number limitation 
and accurate space resolution by the bi-cubic profile. The method is compared with the fourth 
order time-space Yee's scheme in terms of numerical accuracy at nodes and required operation 
counts given accuracy. The numerical convergence rate test shows the method is nearly fourth 
order and the operation counts are comparable with those for the fourth order Yee's scheme 
in conventional computing. With distributed computing implementation the method becomes 
much efficient. Also, the method has the build-in bi-cubic interpolation and thus provides sub- 
grid resolutions at each cell. 

If we use the method with CFL= 1, as shown in Section [4] the method preserves sharp 
profiles in the solution very accurately without any smearing and distortion. The stability and 
numerical accuracy are analyzed. Also, the method computes directly the physical quantities 
e.g., current and electric field gradient, very accurately. The building block of our method is 
the exact integration rule for the Poisson formula against the polynomials. Thus, one can use 
various polynomial approximations locally, including the one using the solution values only. In 
the paper we only implement the methods for the periodic boundary condition but it can be 
extended to the various boundary conditions and class of absorbing boundary conditions. 

Also, we can also extend the interface treatment in [11] for piecewise cubic interpolation at 
the material discontinuity and develop the immersed interface method for discontinuous media. 

An outline of our presentation is: in Section 2 a CIP scheme is proposed, in Section 3 the 
stability and error analysis is presented. Finally in Section 4 we present our numerical tests and 
numerical convergent rate. 
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2 Derivation of the multi-moment scheme 

Consider the two dimensional (TE) Maxwell's equation for magnetic field H = (0,0, H^) and 
electric field E = {E^, Ey,0): 



dE^ _ IdH^ dEy _ IdH, dH, _ 1 fdE^ dE, 



dt £: dy ' dt e dx ^ dt /j. \ dy dx 



(1) 



where the material coefficients e, /j are constants. Let c = be the speed of light. Our 
method uses the fact that ([T]) is equivalent to second order wave equations: 

dfH - c^AH, = 0, d'^E^ - c^AE^ = 0, d'^Ey - c^AEy = 0. (2) 

under the assumption that div(Ex, Ey) = 0. Similarly, the two dimensional acoustic wave 
equation for pressure p and velocity v = {vx,Vy): 

dvx _ dp dvy dp dp _ f dvx , dvy 



dt dx ' dt dy^ dt \dx dy J ^ 
can be treated in exactly the same manner. 

2.1 Poisson's formula and the multi-moment scheme 

Let u{t, xi,X2) be a solution of the wave equation: 

dtu - c^Au = 0, with u(0) = g, ut{0) = h. 
By the Poisson's formula |5], we have 



lB(x,ct) ((ct)2 - 1^ - x|2)2 



27rct 



Using change of variable — x to ^, the solution u{t + At, x) is 

2TTcAtJB ((cA<)2- |^|2)l 27rc7B ((cAi)2- |^|2)^ 

= L((l + e • V)w(t, X + ■)\B) + AtL{dtu{t, X + ■))\B), 
where B = B{0, cA) and 



2vrcAt Jb ((cAt)2 - |^|2)! 
for a function u{t,S.) and ^ = ('^i,'^2)- 

The derivatives dxu{t, x), a = (ai, 02), ai, 02 G N, also satisfy the wave equation 

d^{d^u)-c^A{dy) = o, 
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as long as the solution is smooth, and thus the higher order derivatives of the solution is advanced 
via 

dy{t + At,x) = L[{1 + ^ ■ V)d^u{t,x + ■)\B] + AtL[d^dtu{t,x + ■)\B]. 

We obtain the exact integration formula for the solutions and its derivatives of Hz{t, x), Ex{t, x) 
and Ey{t, x) at x = 0: For a = (0, 0) and a e N^, 

d^H,{t + At,x) = L[{1 + e • V)d^H,{t,x + ■)\B] + AtL[d^dtH,{t,x + ■)\B] 
= L[(l + C • V)d^H,{t, X + -m + fL[df{d^,E^{t, X + •) - d^,Ey{t, x + ■))\B], 

(3) 

d^E^it + At, x) = + e • V)d^E,{t, X + ■)\B] + fL[d^,d^H,it, x + ■)\B], 

d^Eyit + At,x) = L[{1 + ^ ■ V)d^Eyix + ■)\B] - fL[d^^d^H,{t,x + -)\B]. 

Here we use Q to exchange the time derivative and spatial derivatives. 

Let us define a grid of points in the {t, x) space. Let At and Ax be positive numbers. The 
grid is the set of points {tn,Xij) = {nAt,iA,jA) for arbitrary integers {n,i,j). We let d^Hz^j, 
d"Ex2j and d"Ey"j stand respectively for the approximation to the solution d"Hz{tn,Xi,yj), 
d'^E^{tn,Xi,yj) and d^Ey{tn,Xi,yj) for a E (0,0) UN. 

The basis idea of the multi-moment scheme is to define a higher order polynomials P{x, y) 
on each cell [xi,Xj+i] x [yj,yj^i] using grid values including spatial derivatives at four corners 
of the cell, and substitute them to the exact time integration formula ([s]): We evaluate the 
integrals of the polynomials over the ball -6(0, cAt). Henceforth we assume that cAt < Ax, and 
thus the four polynomials are involved in the integration over the ball, and thus the method 
uses variables at 9 nearest grid points. 

We can derive various multi-moment schemes on the basis of ([s]). The resulted scheme de- 
pends on the number of unknown variables we employ at each grid and the order of interpolation 
polynomials; for instance, if we take the grid values, the firs-order derivatives and their second 
order mixed derivatives as unknown variables, we use the bi-cubic Hermite interpolation, known 
as Boger-Fox-Schmit element in finite element methods, 

3 3 

The coefficients of the polynomial are defined by the usual interpolation condition. The resulted 
scheme is written in terms of the bi-cubic polynomial. Let Hz,ij, ^x^ij and Ey -j denote the 
bi-cubic polynomials defined in the cell [xi,Xj4.i] x [yj,yj^i] by the interpolation condition: 

9xHz,ij{x£m) = dx^z^erm d^Ex^ij{x£rn) = 9" E^'^i^, d'^Ey ^^■{xfrn) = ^^Ey^^^, 

for a G {(0, 0), (1, 0), (0, 1), (1, 1)} and xim £ {xij,Xi-ij,Xi^ij^i,Xij^i}. Let us number four 
cells surrounding a gird Xij counter clockwise; Ci = [xj,Xj+i] x [yj,yj^i], C2 = [xj_i,Xi] x 
[yj,yj+i], C3 = [xi-i,Xi] X [yj_i,yj] and C4 = [xi,Xi+i] x [yj_i,yj]. We also number the 
polynomial defined on each cell accordingly, i.e., H^^ij = Hz^i, etc. We let B}^ stand for 

Bi = B{0, cAt) n {x > 0} n {y > 0}, B2 = B{0, cAt) n {x < 0} n {y > 0}, 
^3 = B{<d, cAt) n {x < 0} n {2/ < 0}, Bi = B{0, cAt) n {x > 0} n {y < 0}. 



4 



With these notation, we obtain the CIP scheme: 

(4) 



for a = (0,0) (function value update), a = (1,0) {xi derivative update), a = (0,1) {x2 
derivative update) and a = (1,1) {xi,X2 second order mixed derivative update). As detailed 
in the following sections, Q develops the moments update at time step i„+i based on the 9 



nearest grid moments at time step tn, i-e., update (10). 

One can reduce the number of unknowns at a grid point; for example, one uses the grid 
values, the firs-order derivatives as unknowns. The number of unknowns to be determined at 
each grid point becomes 9: each component H^, Ex,Ey has 3 unknowns at a grid point. Possible 
choices for the interpolation are 

3 k 3 k 

Ck/x'^'^'y'" + c-iAX^y^ + C4,3X^y^ or Ck/x^~'^y^ + c^^ix^y + ci^xy^. 

fc=0 1=0 k=0 1=0 

The coefficients c^/ are determined by using the grid values, the firs-order derivatives at the 
four corners of a cell. The second mixed derivatives being not used, the resulted schemes have 
less complexity than the bi-cubic Hermite polynomial based scheme, however, they produce less 
accurate numerical solution, and the CFL number is less than 1. As for the other choice, we 
consider the bi-liner interpolation: 

co,o + cifiX + co,iy + ci,ixy. 

We then obtain a derivative free nine point scheme. 

2.2 Bi-cubic interpolation and the integration 

In this section, we examine the details for computing the integrals in Q when the bi-cubic 
interpolation is used for the interpolation, i.e., the number of unknowns for Hz, and Ey is 4 
respectively at a grid point; the grid value, the firs-order derivatives and the second order mixed 
derivative. 

2.2.1 Notation 

Let us introduce some notation. For the numerical quantities fij, dx^fij, /ij and 0x^^x2 f hi 
given at the node Xij, we define a 4 x 1 vector (a multi-moment vector): 

— \_fi,jj9xifij,dx2fi,j'>dx-^^x2fi,j'\ ^ ^ ' 

Let P'^ denote a vector composed of the multi-moment vector assigned at the four corner of the 
ceU Cij = [xi,Xi+i] X [yj,yj+i]: 



fT fT nT nT 
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Denote 



X2-X2A 



construct a bi-cubic polynomial 



where di^i = xi^j+i — xi^i and d2j = X2j+i 



X2 



J- 



We 



k=0 £=0 



Xi - Xi 



d 



X2 - X2 



d 



■2J 



di 



(5) 



on the cell Cjj, where the coefficient vector q^'^ are ordered as 

Q'*'"' := (Q'0,0) 90,1' 90,2, ^ZCS) 91,0; f?!,!) ?1,2, 91,3; 92,0> 92,1 , 92,2 , 92,3 > 93,0) 93,1; 93,2, 93,3)^- 

and e(x) denotes 1 x 16 row vectors 

e{x) = [l,Xi,xl,xf\ (g) [l,X2,X2,X2] , 

i.e., the components of e{x) are 

eLm{xi,X2) = x{x2, 

and are ordered as in 

The coefficient q'^'^ of the bi-cubic polynomial Fij is determined by 16 interpolation conditions 
at four corners Xij, Xj+ij, a;j,j+i and Xj+ij+i of the cell: 



dx^ — dx^ fi,j i ^X2 i-^ij) — ^X2 fi,j , 

Fij{Xi + lj)=fi+lj, dxj^Fij{Xi + ij) = dx-^fi + lJ, 8x2^1,] {=Ci + l,j)=dx2fi + l,i, d^ix2^id(^i+^,i) = 9xiX2-f^ + 1->3' 

Fi,j{xi+i,j+i)=fi+i,j+i, 9j;-^Fi,j(xi+ij+i)=9a;j/i+i,j+i, dx2Fij{xi+ij+i)=dx2fi+i,i+i, dx-iX2^i,i(^i+'i-,i+^)='^x-iX2fi+'i-,i+''-^ 

P'i,j{^i,j + l)=fiJ + l, dxiFij(Xij + i) = dx-^fij + l, dx2Fij{Xij + -i)=dx2fi,j + l, d^iX2^id(^i,i + ''-) = ^x-^X2fi,j + 'i-- 



We obtain then the coefficients g*'-' of Fi 



where Q = [Qi, Q2,Q3, Qa] is the interpolation matrix: 



Qi 



r 1 








1 








1 





-3 





-2 





2 





1 








1 

















1 
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-2 
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-3 


-2 
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-2 
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-6 


-4 


-3 


-2 
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1 
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1 


-6 


-3 


-4 


-2 
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1 . 



,Q2 
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,Q3 = 
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-1 
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-3 


-3 
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-2 
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-1 




-6 
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-1 
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-2 
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-3 
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-2 




-6 


3 


2 


-1 


.-4 
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-2 


1 . 




L 4 


-2 


-2 


1 -1 





3 0-10 

-2010 


3 -1 
0-201 



-9-6 3 2 
6 4-3-2 


6 3-2-1 

-4-2 2 1 J 



And Rij is a tensor product of the 4 by 4 identity matrix / and the diagonal matrix with the 
diagonal entries [1, Axi_j, Ax2j, Axi^iAx2j]- 



Ri 



di 



d2, 



di,id2,j 
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Thus we obtain the bi-cubic polynomiah 



dij 



+ e 



dij 



+ e 



X — X 



d 



Q^Rijfij^i. 



Next, let us introduce some matrices for basic operations. For a cubic polynomial p{x) = eo{x)a 
for X G M, where eo{x) = {l,x,x'^,x^) and a = (ao, ai, 02, «3)^, we have 

—p(x) = eo(x)Da, x—pix) = eo{x)MxDa, piax) = eo(x)Daa, pix — s) = eo(x)Tsa, 
dx dx 



where 



D - 



' 


1 





" 




" 








" 




" 1 

















2 





, M = 


1 











= 





Q 

















3 





1 







































1 
















q3 



1 ~x -x-^ 

1 -2x 3a;^ 

1 -3a; 

. 1 

Below, we will use the commutative properties: 

DDa = aDaD, MDa = a-^D^M, TsD^ = DaTs. (6) 

2.2.2 Computation of L 

Now we express the integral in Q in terms of the grid values. We compute the integrals 

L{AFk{xi,j + -)\Bk), 
for k = 1, 2, 3, 4, where A denotes one of the operators 

for ai = 0, 1, a2 = 0, 1, and we renumber the polynomial; -Fi = Fjj, F2 = -F3 = -Fi-ij-i, 

and F4 = Fij-i. Let us denote the corresponding matrix representation for A by Ta, i.e., 

'^l 5i 52 ^2 5i 52 

For the compact expression, we also number the multi-moment vector accordingly, i.e., = 
= f^~^'^ , = P~^'^~^-, and = /*'-'~^. Then from the representations 



"2,j 



i^3(x^J + = e(6 ® (r_i r_i)Q/?i_ij_ir , 

V ''l.i '^2,3-1 J 
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and hence we have 

L{AFi{x,^, + ■)\Bi) = L{e\Bi)TA (d^^D^) QR,,,f 
L{AF2ix,,, + ■)\B2) = L{e\B2)TA (d ^ 



L{AFi{x,^j + ■)\Bi) = L{e\Bi)TA (d 



D_ 



'^2,j-l 



(r_i®T_i)Qi?,„i,,_i/3 

(/®T_i)Qi?,,,_i/4. 



Thus the computation of the integrals is reduced to the computations of 



L{ei^m\Bi) 



1 f e^,m(0 
27rcAt Jb^ ((cAt)2 - |^|2; 



1 dc c?c^ c^c'^ dc dc'^ d^^ dc^ dj^ dc'^ dj" 2dc^ 



Using change of variable ^ = (^1)^2) = cAir (cos0,sin0), the integrals are evaluated as the 
function of := cAt: 

di,e L(e|Si) 

d2,e He\B2) 

d3,c := L{e\B3) 

d4,c := L{e\B4) 



4' 8 ' 12 ' 16 ' 8 ' 67r ' 32 ' 157r' 12 ' 32 ' 60 ' 96 ' 16 ' IStt' 96 ' IOStt 
1 dc 4^ dc^ dc dc^ dc^ d/ d^^ d^^ d^" dj^ d^^ dc^ d/ 2dc 



4' 8 ' 12 ' 16 ' 8 

1 dc dc^ dc^ 
4'~¥' T2''~l6''^ 



67r ' 32 ' IStt' 12 ' 32 ' 60 ' 96 ' 16 ' IStt' 96 ' 1057r_ 

dj' dc^ dc" dj" 2dc^" 
' 96 '"le"' 15^'~ 96"' IOStt 



d_c d_l _d_l d^ d_l _d_l d_l 
8 ' 67r ' 32 ' IStt' 12 ' 32 ' 60 



1 dc dc^ dc^ dc dc^ dc'^ 



dc*" dc^ 



dc^ dc" 



dc^ dc^ 



dc dc 



2d 



6 n 



8 12 



16 8 67r ' 32 ' IStt' 12 



32 ' 60 



96 16 IStt 96 IOStt 



Using these vectors, the integrations are expressed in terms of the vectors and the matrices, 
i.e., we obtain 



L{AFi{x^^j 
L{AF2ix,^j 



■)\Bi) 
■)\B2) 



d2,c Ta iD^_ ® (T_i ® /)Oi?,_ij/ 



= U3,c 



Ta D 



L{AFi{x,^j + ■)\Bi) = d^^^TA [D^®D^ 



From ([6]) , we have for di , ^2 > 



We compute 
di c <S) 



D 



1 ^2 ^3 ^ A/i A^^ ^^3 ^2 ^2^ ^2^2 ^2^3 ^3 ^3^ ^3^2 



4' 8' 12' 16' 8' 67r' 32 ' IStt' 12' 32 ' 60 ' 96 ' 16' IStt' 96 ' IOStt 
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where A = |^ and ji = We denote the right hand side by A(A, fi), i.e., 

-1 ^l fi^ fi^ X Xfi Xfi^ A2 A2/i AV' A^ A^^^ 2a3^3 



A(A,m) = 
Let Ai = A4 = 
die 



12' 16' 8' 67r' 32 ' IStt' 12' 32 ' 60 ' 96 '16' IStt' 96 ' IOStt 



cAt 
dij' 

D 



A2 — A.-? 



/^i = ^2 = ^ and /i3 = ^4 = af^. Then 



_) =A(Ai,^i), d^.^(D^_®D^) =A(-A2,Ai2), 

1 ) = A(-A3, -Ma), d4.c f-D^ «) D 1 ) = A(A4, -^4) 



Therefore we have for ^ = (1 + ^ • V)d'l^d'^^, 



Y,L{AFk{x,^, + -)\Bk 



A(Ai,^i)^ ^1 , A(A2,/i2)T, ^ n^D f2 



k=l 



2j 



"l,i-l"2j- 



"l,i-l"2 

Let US define 16 by 16 matrices Ai, A2, A^, A4: 



jai JQ2 

"l.i": 



2,i-i 



A(Ai,A*i)T^gi?.j 

TaQRij 
TaQR. 



A(Ai,^i) 
A(Ai,^i) 

'^2,j 

A(Ai,Ati) 



^2 



A(A2,Ai2)Tl4(5^ij 



A, = 



A(A3, ^ia) T^(r_i ® T_i)Qi?,_i,,_i 

TA{T^l®T^l)QR^-l,J-l 



A(A3,M3) 
A(A3,M3) 

''2.3-1 
A(A3,At3) 



^4 = 



L "l,i-l"2,i-l 

and let us define 4 by 4 matrices az\ /c = 1, . . . , 9: 



= 1 : 4, ^2 = 5 : 8, ts = 9 : 12, f4 = 13 : 16, 



k{Xi, ^ii)TA{I ®T^l)QRr,,^l 

^^^^^Ta{I®T^i)QR^,j-i 
^^^^Ta{I®T_^)QR,,,-i 



TA{I®T^l)QRr,j^l 



(7) 



A3 

Ai 
A^ 



,ti), 4^ = A3{:,t2)+A4{:,ti), 4'^' = A4(:, t2), 4'^' = ^3(1, ^4) + ^2(:, ^i), 
, h) + A2{:, t2) + vl3(:, ts) + Ai(:,U), 

,t2) + A^{:,t3), a^^ = A2{:,U), 4^ = Ai{:,U) + A2i:,h), 4'^' = ^3). 



Here we use the Matlab notation to express a sub matrix of a given matrix, for instance, for 4 
by 16 matrix A and the index = 1 : 4, we denote by A{:,ti) the sub matrix [^ij], i = 1, . . . , 4, 
j = 1, . . . , 4 of the matrix A. 

Let us label the multi-moment vectors at nearest nine grid points: 

fl = = fi,j-l, f3 = ft+l,j-l) ^4 = f«-l,i, fs = fj.ji fe = fi+l,i) f? = fj-lj + l; fs = fi,j+li fg = fj+l,i+l! 

Then the left hand side of ([T]) for (01,02) = (0, 0), (1, 0), (0, 1), (1, 1) are expressed using the 
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nine vectors and the matrices a'jl'^: 



^L((l + ^V)Ffe(ari,,• + •)|Bfc) 

k=l 

4 

fc=i 

4 

5]L((l + ^V)%F,(xi,,■ + •)|i?fc) 

fc=l 

4 

J2L{{l+^-V)dl^^F,{x,,, + -)\Bk) 



fe=i 



fe=i 



As for A = d^.d^^dg^ and A = d^^dg'dg\ we obtain 



k=l 



A3A(A3,^3) 



"l,i-l"2,j 



j-1 



"l,i-l"2,j-l 



TA{I®T_^)QRr,,-lf 



and 



cAt^L(AFfe(a;,,j +-)|-Bfe) 



/^iA(Ai,/xi) 

/^3A(A3,/X3) 

jai 702 
"l,i-l"2,j-l 



jai ,a2 
"l,i-l"2,j 



"l,i"2,j-l 



Similarly, one has 

4 

fe=i 

4 

fe=i 

4 

fe=i 

4 



fe=i 



Y^^k^k, cAt 



fe=i 



4 

fc=i 

4 

fe=i 

4 

^L(%a|j,Ffe(ar,,, + -)|Bfe) 



fe=i 



(8) 



fe=i 



(9) 



with 4 by 4 matrices 6^'-', Cj^.'-^, A; = 1, . . . , 9 which are defined in the same manner. 
2.3 The multi-moment scheme 

Suppose the numerical approximations to the exact solutions and its first derivatives 

dx^H'^ixij) and dx2H'^{xij), and the second derivative d'^^^^-^zi^id) known at all grid 
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points Xij = {xi,yj) at time step tn, which we denote by 

d /?" d /?" (9^ ft" 

Similarly, for numerical approximations to Ex and Ey and the derivatives, we use symbols e" j ^ 

and e!^,' Let hf •, exf,- and ey",- denote the multi moment vectors at grid; 

yi''ij ''ij ''tj '■ij 

^hj — l"'zp'Jxiii'tj,Ux2"'ijj'-'xiX2"'ij\ ' 
„n _r„nn„nfl„nfl2 f> n]T 
^xij — y^xij 1 ^xi ^xij , ^X2 ^xij ) ^xiX2 ^^iji ' 

*^yij ~ l^yij^^^i^yiji^x2^yij^^xiX2^yiji ■ 

We number the numerical solutions in the way as above. Based on the exact integration formula 
(|3j), and ([8]),(|9j), we arrive at the multi-moment scheme for the Maxwell's equations: 

n+l _ Jdun , 1 Jd ^ n 1 n 



un+1 _ spy „»Jun, J_«J n _J_i. 



ev 



n+l _ v^9 1 i,j 



ELi^4"K + <"e4, (10) 



n+l _ v^9 _J_J,«Jun -Un*'-?! 



Thus, the method uses 9 nearest neighbor points for 12 components (4 moments) for H^, E^, 
and Ey. Each 4 x 36 matrices [a*^'-', . . . , Og-'], [6*^'-', . . . , 6g-^] and [c^'-', . . . , Cg-'] has 100 nonzero 



entries. Hence the total cost for the update ( 10 ) amounts to 700. 



We would like to emphasis that the multi-moment scheme provides with the first and the 
second derivatives as we as the function value at each grid. When displaying the numerical 
solution, one can construct the bi-cubic polynomial and can evaluate any spacial point with 
using the interpolation. 

3 Stability 

We analyze the stability of the multi-moment scheme. For simplicity, we assume that the grid 
length is uniform, i.e.. 

Ax := Xi - Xi-i = yj - yj-i 



for all In this case, the 4 by 4 matrices a!'-', and c!'"' in (10) remain the same throughout 



of (^,j), and thus we omit the subscript. Suppose that f = {fe,j}e,j is N periodic with respect 
to i.e., 

fo,j = ^NJ, f£,o = ^e,N, 
for all < I, j < M. Let us consider the discrete Fourier transform: 

N N N N 

^ ^ ^ij^ ~''^^''^^[fij^'^xifi,j,dx2fi,j,dx-^x2fhj] ^ 

i=0 j=0 i=0 j=0 

where k ■ Xij = {ki,k2) • {xi,yj) = kiXi + k2yj- The discrete Fourier transform of the sequence 
Fij := aifi_ij_i+a2fjj-i+a3fj+ij_i+a4fi_ij+a5fij+a6fi+ij+a7fi_ij+i+a8fjj4_i+a9fi+ 
becomes: 

Af N 

J^F„e-2-^^-"^.i =5({afe},2^A:iAx,27rA:2Ax)f[A:], (11) 

t=0 j=0 
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where g{{ak},27rkiAx,2'iTk2Ax) is 4 x 4 matrix depending on A, kiAx and k2Ax: 



g{{ak},27rkiAx,27:k2Ax) = e^^'^-^^-^^^^^'m + e^^^^^^^Aa;^^ ^ ^2^i(k^-k2)^x^^ 
+ e-2^^^i^^a4 + as + e^^^^^^^ae + e2'^*(-'=i+'=2)^^a7 + e^^^^^^Ax^^ ^ g2^i(fci+fc2)Ax^g^ 

Now let us assume {Hz, E^, Ey) is periodic. Let h" = {h"j}jj, ex" = {ex"j}jj and ey" = 
{ey"^.}jj. From (10) and ( [Il| ), the amphfication factor <I>(A, 27rA;iAx, 27rA;2Ax) (12 by 12 matrix) 
of the proposed scheme is given by 



H\0i,e2) 



g{{ak}, 01,62) 

^9{{ck},di,d2) 9{{ak} 



^^gi{ck},eu02) 



3 

ce 



n,v2) 



o 

g(.{ak},0: 



2) 



(12) 



where 61 = 27rA;iAx, 62 = 27r/c2Ax, i.e., <^(A, 27rA;iAx, 27rA;2Ax) maps (h",ex",ey") to the next 
step (h"-+i,e. 



n+l 



, Cy 



n+l} 



Let p{X, 61, 62) denote the set of eigenvalues of <1>(A, Oi, 62) ( 9 eigenvalues). Figure [T] shows the 
maximum absolute value of the eigenvalues, max{|p(A, 6*1, ^2)!}, against [6*1, 6*2] G [— vr, 0] x [— vr, 0] 
for A = 0.2, A = 0.5 and A = 1, respectively. Numerically we find that all eigenvalues have the 
magnitude equal to or less than 1 for arbitrary {61,62)- The magnitude is close to 1 in a wide 
range of [^1,^2], which indicates that the numerical scheme is less dissipative. 



4 Numerical test 

In this section, we show the numerical performance of the multi-moment scheme through some 

numerical tests. 

Example 1. Plane waves. 

In this example, we compute the numerical solutions for plane waves. We compare the numerical 
solutions with those produced by the fourth order in time and space FDTD (Yee's scheme). The 
Yee's scheme computes E field and H field at different time level, and thus one must provide 
the exact initial condition at time t = for E field and t = for H field to obtain an 

accurate numerical solution. The plane wave solution is suitable to avoid the issue with the 
initial condition since the exact solution is easily obtained. 

2 

Let fa{x,y) = exp(— ^) for o" > 0. Let fa also denote its periodic extension to x direction 
with periodicity L, i.e., /^(x + L,y) = fa{x,y) for all {x,y) G M^. We rotate the function with 
the angle 6 with respect to the origin (0, 0) to construct the one way propagating plane wave 
solution for Maxwell's equation: 

Hz{t,x,y) = fa{xcos6-ysm6-t), E^{t,x,y) = Hz{t,x,y)sm6, Ey{t,x,y) = Hz{t,x,y) cos 6 

If we set L = cos 6 with 6 = tan^^ m for m £ N, they are the periodic solution of the Maxwell's 
equation for e = /i = 1 in the domain Q = [— |, 5] x 5]- In this numerical test we consider 
the solutions 

3 

{t,x,y) = ^ fa {x cos 0,ri - y sin 6,n - t) , 

Tn—0 

3 3 

Ex{t,x,y) = ^ H^{t,x,y)sin9m, Ey{t,x,y) = ^ Hz{t,x,y) cos9,n, 

m— m— 
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where 6m = tan~^ m. We test the case a~'^ = 500, 1500. Figure 1 shows the initial profile of 
Hz{0, x,y), and the four plane waves at t = 0, fa{xcos6m — ysinOm) for m = 0,1,2,3, with 
0"^^ = 500, in the domain. The arrow in each plot shows the direction of wave propagation. 

We report the accuracy of the multi moment scheme. The time step size is fixed to be 
At = AAx for each mesh size Ax = A^~^, G {50,100,150,200}. We report the numerical 
solutions for A = 1 and ^ = The numerical solutions at time T = 1 are produced by the 
multi moment scheme and compared to the exact solution. The number of iteration is N for 
A = 1, and IAN for X = where the numerical solution approximates the solution at time 

T = ^ ~ 0.98995. 

The initial value for h^j is provided by the exact solution: 

hlj ^ H^{0,Xi,yj), da:h°.j ^ d^H^iO, x^y^), dyh°.j ^ dyH^{0,x^,yj), d^yhl^ = dlyH^{0,x„yj). 

Similarly, ^ and j are given by the exact solution. In a practical situation, the initial 
condition in function form may not be available, and we only have the function value at each 
grid point. In that case we employ finite difference of the initial grid function to provide the 
initial condition for derivatives. 

For each mesh size j^, the error in the numerical solutions is measured by norm: 

ei = m^x{\hl^ - H,{T,x.,,,)\, |(e,)^, - E.,{T, x,.,)\, |(e,)I^ - Ey{T,x,,,)\} 

The multi-moment scheme computes the first derivatives and the second mix derivative as 
numerical solutions. We report the relative error in the first derivatives: 

_^ \do,hl^~-d^H,{T,x,,j)\ \{d^e.,yi^-d^E,{T,x,,,)\ \{d^ey)lj - d^Ey{T,x,^,)\ 

\do.H,{T,x,^,)\ ' \dMT,x,.j)\ ' \do.Ey{T,x,,j)\ ^' 

where a = x,y. 

For comparison, we also report the numerical solution produced by fourth-order in time and 
space FDTD (see [1]). In the FDTD numerical simulation, we employ the CFL number to be 
with which the FDTD method provides the best performance, i.e., FDTD with CFL= 

produces most accurate numerical solutions among the other CFL. In this test, the initial values 

_ 1 _ 1 

for h^j, exj^j and ey^^ are given exactly: 

_ 1 _ 1 

= H^{0,Xi,yj), {ex)ij = Ex{-^,Xi,yj), {ey)^J = Ey{-^,Xi,yj). 

When the analytic solutions are not available, one must solve the Maxwell's equations backward 
in time to provide an accurate initial condition for Ex{—^,Xi,yj) and Ey{—^,Xi,yj) to start 
the FDTD scheme, by employing another time marching method such as Runge-Kutta schemes. 
We compute the error in the numerical solution: 

ei = max{\hl^ - H,{T,x.,,,)\, |(e,)^^ - E.,{T, x,^,)\, |(e,)^, - Ey{T,x,.,)\}, 

\dahlj-dMT,x,,j)\ \{d^e^)lj-dMT,x^,j)\ \{dc,ey)lj - do^EyjT, x,^j)\ 

i,j.aex,y \daH^{T,X,j)\ ' \daEx{T, Xij)\ ' \daEy{T,Xij)\ 

where dah^j, {daex)^j and {daey)^j are computed by employing third order finite difference 
which is used in the fourth order FDTD scheme. 
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Figure [3] shows the error ei and £2 of the solutions generated by the multi moment method 
with A = 1 and = ^-^d the fourth order FDTD with X = for the initial profile with 

o"^^ = 500 against the grid number N. The order of accuracy of the numerical solutions for each 
method is shown in Table [l] And Figure |4] and Table [2] are numerical results when cr"^ = 1500. 

Let us consider the total computational cost in the multi-moment scheme to obtain the 
numerical solution at T = 1. For each time step, 300 A^^ operation is required to update Hz, and 
200 for gach Er, and Ey. Since the number of time integration in the multi-moment scheme 
with A = 1 is the total cost (operation count) amounts to 700 A^^ x N = 700N^. 

Let us focus on the numerical solutions by the multi-moment scheme and Yee's scheme when 

A^ = 50. Figure [3] shows that the error in the numerical solution produced by the multi- moment 

scheme with A = 1 is 10~^ while the error by FDTD is 10~^. So, the multi-moment scheme is 

10 times more accurate than the FDTD with A = ^ for this mesh size. To obtain the same 

V 2 

accurate numerical solution by FDTD, we must take the half mesh size 1/(2A^) = 1/100. For 
the fourth order Yee's scheme, the cost for the one step map is 51 and the number of time 
integration is ~ 2^/2N, thus the total cost for the numerical solution at T = 1 amounts to 
2V2N X 51(2Ar)2 ^ 577^r3. 
Example 2. Sharp profile. 

Next we solve ([T]) with e = /x = 1. The initial condition is 



H{0,x,y) = \ J' , E^{0,x,y) = Ey{0,x,y) = 0, 



where D = [0.25, 0.75] x [0.25, 0.75]. The mesh size is 0.01. The initial condition is approximated 
by bi-cubic polynomial with the first and second derivatives being 0. We do not use any other 



techniques to approximate the initial discontinuous profile. We use our algorithm (10) with 
CFL= 1. The initial condition and the numerical solution for H at time T = 0.15 and T = 0.25 
are displayed in Figure [5j No oscillation is observed in the numerical solution. We also employed 
the fourth order FDTD with the same initial condition. Numerical oscillations were found near 
the sharp profile. 

Example 3. Hidden resolution 

We solve ([T]) with e = /i = 1 with the initial condition is 

H{0,x,y) = exp (^-^^^ ' Ex{0,x,y) = Ey{0,x,y) = 0. 

The mesh size is ^ and CFL= 1, and so At = ^. We compute the numerical solutions 
at t = WAt. We denote the numerical solutions by {hl^ {9xhl^ i j} i'^y^l^ i 

{d^yhl^.^}ij. When visualizing the numerical solutions, we usually construct the bi-linear 
interpolation in each cell using the numerical solution at the grid. For instance, for the vi- 
sualization of the solution {hl^ we plot the bi-linear interpolation constructed by using 
the gird value {hl^ and for visualizing the numerical solution {dxhl^ ^ we plot the 
bi-linear interpolation constructed by the gird value {dxh]P These two bi-linear interpola- 
tions are unrelated. In the left column in Figure [6j we plot the piece wise bi-linear interpolation 
for {hl^^j}ij (top) and the one for {dxhl^^ (bottom). One can observe the spiky peaks and 
dips in the plots. 

As have been mentioned repeatedly, the multi-moment scheme produces the derivatives as 
well as the function value at each grid. This is equivalent to state that the multi-moment scheme 
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computes the bi-cubic polynomial in each cell as a numerical solution. So when plotting the 
numerical solution, we should use the computed bi-cubic interpolation instead of bi-linear inter- 
polation. Let us construct the bi-cubic polynomial hz^ij{x,y) in each cell x [yj-i,yj] 
using the numerical solutions, and let hz{x,y) denote the piece wise bi-cubic polynomial defined 
in the domain [—0.5, —0.5] x [-0.5, —0.5], i.e.. 



hz{x,y) = hz,ij{x,y), {x,y) £ [xi-i,Xi] x [yj^i,yj]. 



(13) 



In the right column of Figure [6| we show the surface plot of the bi-cubic interpolations hz{x, y) 
(top) and dxhz{x,y) (bottom). The numerical solutions are depicted as smooth functions. 
Derivative free method 

We have also implemented the method using the bi-linear interpolation at each cell Cij: 

F{x, y) = F(0, 0)(1 -x){l-y) + F{Ax, 0)x(l - y) + F{0, Ax){l - x)y + F(Ax, Ax)xy 
In this way we obtain a derivative free nine point scheme: 



— L/o — — L'i 
Li 
Li 



where 

L2 



A2 

A _2-A 
^8' 4^ 



2A)A A2 (7r-2A)A 
' 2^' "hi 



A A 2- 



1 - 2A 

A A 

^ 8 



2A2 (7r-2A)A A^ (tt - 2A)A A^ 
~' 27r ' 27r' 2tt ' 2tt 
2-A 2-A 



A A 



A „ A 



This method is also stable with CFL = 1 but is second order accurate. Our numerical tests 
show that if we let A = 1, (CFL= 1) then there is no significant dissipation but A = 0.5 it has 
30% dissipation at T = 1 with speed one. For oblique plane waves there is no significant phase 
error with CFL = 1. An advantage of this method is that it is simple to be implemented and 
to be extended to the three dimension case. 



5 Conclusion 

We developed a numerical method for solving time-domain Maxwell's equation. It is fully explicit 
space and time integration method with higher order accuracy and CFL number being one. The 
bi-cubic interpolation is used for the solution profile to attain the resolution. It preserves sharp 
profiles very accurately without any smearing and distortion due to the exact time integration 
and high resolution approximation. The stability of the method were analyzed, and the nearly 
forth order accuracy were observed. 
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Figure 1: The maximum absolute value of the eigenvalues of the amplification matrix (12) 
against (6*1,612). A = 0.2 (left), A = 0.5 (center), A = 1, (right). 




Figure 2: Initial profiles Hz{0,x,y) = Ylm=o f<^(^^^^^'"T- ~ 2/ sin 0m) (left) and fa{x cos9m — 
y sin 9m) for m = 0, 1, 2, 3. The inverse of the variance a~'^ is 500. The arrow in each plot shows 
the direction of wave propagation. 



Table 1: Numerical errors ei and 62 and the order of convergence, a in the initial profile is 



a-2 = 


500. 

Multi moment (A = 


= 1) 


Multi moment (A = 


A) 


FDTD (A = ^) 


N 


50 


100 


150 


200 


50 


100 


150 


200 


50 


100 


150 


200 


ei 


9.86e-3 


8.93e-4 


1.88e-4 


6.06e-5 


1.07e-l 


1.72e-2 


5.45e-3 


2.32e-3 


1.04e-l 


7.15e-3 


1.37e-3 


4.39e-4 


order 




3.46 


3.84 


3.93 




2.63 


2.84 


2.97 




3.86 


4.06 


3.96 




1.72e-2 


1.50e-3 


3.20e-4 


1.04e-4 


1.30e-l 


2.20e-2 


6.86e-3 


2.95e-3 


1.29e-l 


1.07e-2 


2.32e-3 


7.33e-4 


order 




3.51 


3.81 


3.90 




2.56 


2.87 


2.92 




3.58 


3.78 


4.00 
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Figure 3: Numerical error ei (left) and £2 (right) in the numerical solutions for the initial 
condition o""^ = 500. The order of accuracy is shown in Table [l] 




N N 



Figure 4: Numerical error ei (left) and £2 (right) in the numerical solutions for the initial 
condition cr^^ = 1500. The order of accuracy is shown in Table [2] 



Table 2: Numerical errors ei and €2 and the order of convergence, a in the initial profile is 
0--2 = 1500. 





Multi moment (A = 


= 1) 


Multi moment (A — 




FDTD (A = ^) 


N 


50 


100 


150 


200 


50 


100 


150 


200 


50 


100 


150 


200 


ei 


8.80e-2 


1.08e-2 


2.70e-3 


9.10e-4 


3.65e-l 


1.12e-l 


4.04e-2 


1.84e-2 


5.23e-l 


9.20e-2 


1.96e-2 


6.46e-3 


order 




3.02 


3.41 


3.79 




1.69 


2.52 


2.72 




2.50 


3.80 


3.86 


£2 


1.63e-l 


2.21e-2 


5.35e-3 


1.89e-3 


6.33e-l 


1.51e-l 


5.54e-2 


2.47e-2 


5.19e-l 


1.43e-l 


3.50e-2 


1.09e-2 


order 




2.87 


3.50 


3.60 




2.06 


2.47 


2.80 




1.85 


3.47 


4.03 
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Figure 5: The initial condition for H (the first column) and the numerical solution for H at 
time T = 0.15 (the second column) and T = 0.25 (the third column). 
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-0.5 -0.5 



-0.5 -0.5 



Figure 6: Left column: 2D surface plot of {hl^ ^ Ai j{top) and {S^/i^'^- jij (bottom). Right 



column: The bi-cubic interpolation (13) hz{x,y) (top) and dxhz{x,y) (bottom). 
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